| A generic implementation of Branch & Bound

> {-# LANGUAGE BangPatterns, ViewPatterns, GeneralizedNewtypeDeriving,
>         MultiParamTypeClasses, FunctionalDependencies,
>         UndecidableInstances, FlexibleInstances, StandaloneDeriving #-}

> module BB (BBable(..), Box(..), Bx, inf, infIO, infU, box2bb) where

> import qualified Data.List as List
> import qualified Bs
> import qualified Aux
> import qualified Maybe
> import qualified Iv 
> import Control.Monad (when)
> import Control.Monad.Trans
> import Control.Monad.Identity

TODO make this type abstract

> newtype Bx a = Bx (a, [(Bool,Bool)])

> newtype B b = B b deriving (BBable a)

> class (Ord a, Num a) => BBable a b | b -> a where
>   branch :: b -> [b]
>   bounds :: b -> Iv.Iv a
>   -- ^ Calls to bounds must be cached.
>   info   :: b -> Int -> String

> lb :: BBable a b => b -> a
> lb = Iv.lo . bounds

> ub :: BBable a b => b -> a
> ub = Iv.hi . bounds

> unB :: B t -> t
> unB (B b) = b

> instance BBable a b => Eq (B b) where
>   B b1 == B b2 = lb b1 == lb b2

> instance BBable a b => Ord (B b) where
>   compare (B b1) (B b2) = compare (lb b1) (lb b2)


TODO handle verbosity through MonadPrinter

| Returns an interval containing the infimum and a minimizer.

> infU :: (Aux.MonadPrinter p, BBable a b) => b -> a -> a -> Int -> p (Iv.Iv a, b)
> infU b0 u eps v = inf_ [B b0] u
>   where
>     inf_ (b:bs) u =
>       do when (v >= 1) $ 
>            Aux.printLn (show (length bs + 1) ++ "\t" ++ 
>                      show u ++ "\t" ++ 
>                      show (lb b) ++ "\t" ++ 
>                      info b v)
>          if u - lb b < eps then return $ (Iv.iv (lb b) u, unB b) else
>            let b' = List.sort $ branch b in
>            let u' = foldr min u $ map ub b' in
>               inf_ (filter ((u' >=) . lb) $ Aux.merge b' bs) u'
>     inf_ []     _ = error "inf error: the problem evaporated"
>                     -- This may happen with floating point numbers on Intel.

| Infimum with known upper bound

> inf_ :: (Aux.MonadPrinter p, BBable a b) => b -> a -> Int -> p (Iv.Iv a, b)
> inf_ b = infU b (ub b)

> inf :: BBable a b => b -> a -> (Iv.Iv a, b)
> inf b eps = runIdentity (inf_ b eps 0)

> infIO :: BBable a b => b -> a -> Int -> IO (Iv.Iv a, b)
> infIO = inf_
  
> class (Ord a, Num a) => Box a b | b -> a where
>   dim     :: b -> Int
>   boxBounds :: b -> Iv.Iv a
>   boxInfo :: b -> Int -> String
>   sharp   :: b -> Int -> Maybe Aux.Hand
>   proj    :: b -> Int -> Aux.Hand -> b
>   pickDim :: b -> Int
>   subdiv  :: b -> Int -> (b,b)

| Try all dimensions beginning at k, stop at first success, otherwise
fail.

> findSharp :: Box a b => b -> Int -> Maybe (Int, Aux.Hand)
> findSharp b k = 
>   let n = dim b in fmap Maybe.fromJust $ List.find (/= Nothing) $
>   map ((\i -> fmap ((,) i) $ sharp b i) . (`mod` n)) [k..k+n]

| Do all possible projections, beginning at dimension k. If
elimination is possible, return Nothing.

<TODO> Give back sequence of boxes, so we can stop as soon as we have
enough information, i.e.:

- the answer is epsilon-tight, or

- are better than the next-best box (this might be an interesting
  application of lazy evaluation)
</TODO>

> deflate :: Box a b => Bx b -> Int -> Maybe (Bx b)
> deflate (Bx(b,p)) k = 
>   maybe (Just $ Bx(b,p)) (\(i,h) ->
>     let (p0,(l,r):p') = List.splitAt i p in
>     if (h == Aux.L && not l) ||
>        (h == Aux.R && not r) then Nothing
>     else 
>       let b' = proj b i h in
>       let n' = dim b' in 
>       if n' == 0 then Just $ Bx(b,p) else
>       deflate (Bx (b', p0++p')) (k `mod` n'))
>   $ findSharp b k

> instance Box a b => BBable a (Bx b)  where
>   branch (Bx (b,p)) = 
>     let k = pickDim b in
>     let (b1,b2) = subdiv b k in
>     let (p0,(l,r):p') = List.splitAt k p in
>     Maybe.catMaybes [ deflate (Bx (b1, p0++(l,False):p')) k
>                     , deflate (Bx (b2, p0++(False,r):p')) k ]
>   bounds (Bx (b,_)) = boxBounds b
>   info (Bx (b,_)) = boxInfo b

> box2bb :: Box a b => b -> Bx b
> box2bb b = Maybe.fromJust $ 
>            deflate (Bx (b, replicate (dim b) (True,True))) 0
